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Crossover phenomenon in self-organized critical sandpile models 
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We consider a stochastic sandpile where the sand-grains of unstable sites are randomly distributed 
to the nearest neighbors. Increasing the value of the threshold condition the stochastic character of 
the distribution is lost and a crossover to the scaling behavior of a different sandpile model takes 
place where the sand-grains are equally transferred to the nearest neighbors. The crossover behavior 
is numerically analyzed in detail, especially we consider the exponents which determine the scaling 
behavior. 
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I. INTRODUCTION 

Crossovers between different universality classes are 
well known from equilibrium phase transitions. Sim- 
ilar crossover phenomena are also known from non- 
equilibrium systems which exhibit self-organized critical- 
ity Q . For instance the forest-fire model of Drossel and 
Schwabl [|| displays a crossover to a percolation like scal- 
ing behavior if one introduces an immunity parameter 
which prevents trees from burning j^] . Another example 
is the directed Abelian sandpile model of Dhar and Ra- 
maswamy B. Introducing a second stochastic toppling 
condition the system changes to a directed percolation 
like scaling behavior || . 

In this paper we investigate a crossover between the 
stochastic sandpile model introduced by Manna H and 
the Zhang sandpile model 0. Following Ben-Hur and 
Biham j| the crossover connects the different universal- 
ity classes of undirected sandpile models (the energy or 
sand-grain transfer of the nearest neighbors is isotropic, 
e.g. the Zhang model) and undirected on average sandpile 
models (the energy transfer is isotropic on average only, 
e.g. the Manna model). Considering the Manna model 
the crossover takes place just by increasing the threshold 
value which determines the dynamics, i.e., no additional 
parameter has to be introduced. 

In the next section we briefly remind the distinctive 
characteristics of the Zhang model. These characteris- 
tics will allow us in the following to identify the typical 
Zhang scaling behavior. Then we describe in section III 
the Manna model for different values of the threshold 
condition. The crossover behavior between both mod- 
els is investigated in section IV. A summary closes the 
paper. 



II. THE ZHANG MODEL 

Consider the Zhang model [jrj on a two dimensional 
square lattice of linear size L. A continuous value E r > 
representing the energy is associated to each lattice site 
r. A configuration {-Ey} is stable if E L < E c for all 
lattice sites r. For the sake of simplicity we choose in all 
simulations E c — 1. A quantum of energy 5 is added to 
a randomly chosen lattice site r, i.e., 



E r 



E.. 



(1) 



We consider in this work especially the slow driving limit 
8 <C E c , For 5 — > all lattice sites grow parallel and the 
Zhang model corresponds to the conservative limit of the 
"spring block" model of Christensen and Olami || . 

In the case that due to the perturbation a site exceeds 
the critical value E c an activation event will occur, the 
unstable site relaxes to zero, and the energy is added to 
the nearest neighbors, i.e., 



E r 



0. 



E, 



r.NN 



E 



r.NN 



(2) 



(3) 



The transferred energy may activate the neighboring sites 
and thus an avalanche of relaxation events may take 
place. Energy may leave the system only at the bound- 
ary. Since the energy of unstable sites is equally trans- 
fered to the nearest neighbors it was argued that the 
Zhang model belongs to the universality class of undi- 
rected sandpile models [||. It was expected that the 
Zhang model and the well-known Bak-Tang-Wiesenfeld 
(BTW) model § belong to the same universality class. 
But the scaling behavior of the BTW avalanches is com- 
plex and is not understood. Although most authors agree 
upon a breakdown of simple scaling the interpretation of 
the numerically obtained data is still controversial among 
the different groups (see for instance |ic|-[l^|). However, 
we use in the following the classification ansatz of j|] and 
denote the universality class of the Zhang model as the 
class of undirected sandpile models. 

The avalanches are characterized by several physical 
properties like the size s (number of relaxation events) , 
the area a (number of distinct toppled sites), the time 
t (number of parallel updates until the configuration is 
stable), the radius r (radius of gyration), etc. In the 
critical steady state the corresponding probability distri- 
butions should obey power-law behavior (T| 



Px(x) 



with 



(4) 



X G 



characterized by the avalanche exponents t x 
{s,a,t,r}. 

The Zhang model was intensively investigated in the 
last years (see for instance j7],[l4] Q). A characteristic 



accepted for publication in Physical Review E 62, (2000) 



2 



property of the Zhang model is the concentration of the 
steady state energy distribution p(E) around z distinct 
peaks, where z is the number of nearest neighbors of the 
lattice Jfl, |l3Jl5| , [[7j . The peaks are located at multiples 
of (z + 1)/? and the height of the peaks diverges in the 
thermodynamic limit L —> oo (see (ItJ and references 
therein). In the case of an infinite lattice the energy 
distribution is given by 



p(E) = £ /< S(E - Ek 



(5) 



i=0 



where fi denotes the statistical weight and E\ denotes 
the position of each peak. It was found numerically that 
the statistical weights are independent of the input en- 
ergy 5 [0. Thus the statistical weights can be regarded 
as another fingerprint of the Zhang model. 

Analyzing the numerically obtained avalanche distri- 
butions [Eq. (|])] it was observed that the avalanche ex- 
ponents of the Zhang model exhibit finite-size corrections 
according to @p| 



t x (L) 



const,. 



(6) 



In this case the values of the infinite lattice r x are ob- 
tained by an extrapolation to L — > oo. More than the 
explicit values of the exponents this characteristic system 
size dependence allows us in the following to identify the 
Zhang like scaling behavior. 



III. THE MANNA MODEL 

A stochastic sandpile model in which integer values 
represent local energies (or number of sand-grains) was 
introduced by Manna ||. Here, unstable sites relax to 
zero if E r > E c and the removed energy is randomly 
distributed to the nearest neighbors in the way that one 
chooses randomly for each energy unit (one sand-grain) 
one neighbor. Thus the Manna model is characterized 
by a stochastic energy transfer and according to || it 
belongs to the so-called universality class of undirected 
on average sandpile models. 

Due to the reduction of the energy of unstable sites 
to zero both the Zhang and the Manna model are non- 
Abelian models jl{|, i.e., the stable energy configurations 
depend on the sequence in which unstable sites are top- 
pled. Recently Dhar introduced an Abelian version of the 
two-dimensional Manna model where the energy of crit- 
ical sites is not reduced to zero but Eij — > E^j — 2. The 
energy AE = - 2 is then equally distributed with proba- 
bility 1/2 to the sites (i ± 1, j) or otherwise to the sites 
± 1) ]20|] . In this case it is possible to extend an op- 
erator algebra, which was successfully applied in study- 
ing the Bak-Tang-Wiesenfeld model , to this modified 
Manna model. 

In contrast to this analytically tractable Abelian 
Manna model we consider in this paper the original 
Manna model (E L — > 0). Numerous numerical analy- 
sis of the Manna model were performed for E c = 2 and 



the values of the exponents are known within some error- 
bars (see for instance ||||,[l^,|lj and references therein). 
Usually one assumes that the value of the critical en- 
ergy E c has no influence on the scaling behavior of the 
model (see for instance H^l)) and indeed the analysis 
of the avalanche probability distributions [Eq. (f|)] for 
E c G {2, 3, 5, 10} reveals that the exponents are inde- 
pendent of E c (see Fig. 1 in ref. ^J). But as we will 
see a crossover to a different universality takes place for 
sufficiently large values of the critical energy. 

Consider the Manna model on a lattice with z nearest 
neighbors. Assume that a given site r exceeds the critical 
value, i.e., E L > E c . Thus a toppling process takes place 
and the unstable energy E r is randomly distributed to 
the nearest neighbors. A given neighbor of r get n energy 
units and the corresponding probability distribution is 



P(n,E L ) = 



E, 



p q 



(7) 



with p = 1/z and q = 1 — p. The corresponding expec- 
tation value of the energy transfer is p = pE r and its 
variance a = ^JpqE r . One can therefore distinguish be- 
tween two regimes: for small values of £V the expectation 
value is of the same order as the variance (p « er), i.e., the 
energy transfer to the nearest neighbors is characterized 
by strong fluctuations. With increasing unstable ener- 
gies the fluctuations decreases and they can be neglected 
if n ^> a for large values of E T . 

Since the value of the unstable energy E L is of the 
same order as the critical energy E c we introduce the 
crossover parameter K = y/Jz — 1) / E c . As long as the 
fluctuations are relevant K » 1 we expect that the scal- 
ing behavior agrees with that of the Manna model for 
small E c . Increasing the threshold value E c the fluc- 
tuations become irrelevant if K -C 1. In this case the 
energy is nearly equally distributed to the nearest neigh- 
bors and the model corresponds to the Zhang model in 
the limit 5 <C E c . Thus we expect that a crossover from 
the Manna to the Zhang scaling behavior takes place if 
one increases the value of the critical energy. The details 
of this crossover are investigated in the next section. 



IV. CROSSOVER BEHAVIOR 
A. Energy distribution and energy correlations 

In this subsection we investigate the static properties 
of the Manna model for various values of the critical en- 
ergy E c . We show that the crossover affects the distri- 
bution of the energies, its scaling behavior as well as the 
spatial correlations of the energies. 

We measured the average energy in the steady state 



(Eh 



(8) 



for various values of L and E c . Figure [TJ shows the sys- 
tem size dependence of the average energy (E)l for three 
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FIG. 1. The average energy in the steady state {E)l as a 
function of the inverse system size L for different values of E c . 
The extrapolation to the vertical axis yields the value of the 
average energy for L — > oo. 



different values of E c . For all considered values of E c we 
observed that the system size dependence of the average 
energy is given by 



(E)l = (E)! 



const/L. 



(9) 



This behavior was already found by Manna for E c = 2 
H and it is also known from the Bak-Tang-Wiesenfeld 
model |^2f . The origin of this system size dependence are 
boundary effects. The average energy on the boundary 
is smaller than the energy in the bulk and the relative 
number of boundary sites scales as L _1 . 

The extrapolation to L — > oo yields the average en- 
ergy of the infinite system size and the obtained values 
are shown in Fig. | as a function of E c . Increasing the 
critical energy from E c — 2 the average energy (E) l^oc 
decreases and reaches a minimum for E c w 100. Here the 
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FIG. 2. The average energy in the steady state (E) l^oo as 
a function of the critical energy E c . The inset displays the 
average energy vs. 1/E C for E c > 350. The corresponding 
value of the pure Zhang model is marked as E c — oo and as 
an open circle (inset), respectively. 



FIG. 3. The energy probability distribution p(E) for 
L = 128 and various values of E c . To guide the eye solid 
lines are plotted instead of symbols. Normalizing p(E) to E c 
allows to illustrate how the peak structure of the distribution 
appears with increasing critical energy. This is shown in the 
lower right figure for E c = 50, 500, 10000. The vertical solid 
lines correspond to the position of the peaks of the Zhang 
model Ei = 5/16 i with i = 0, 1,2,3 @. 



behavior changes and the average energy increases with 
the critical energy. For large values of E c the average 
energy saturates in the vicinity of the value of the Zhang 
model. A detail analysis suggests that the dependence of 
average energy on E c is given by 



(E)i 



(E) zhang - const/£; c 



(10) 



for E c > 350 (see inset of Fig. g). Since finite values 
of the critical energy results in finite fluctuations we get 
from Eq. (|l^) that the average energy is affected by these 
fluctuations even for very large values of E c . 

After the average value of the energy we consider now 
the steady state distribution p(E) of the energies. In 
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FIG. 4. Analogous to Fig. |^ but for a honeycomb lattice. 
The vertical solid lines correspond to the position of the peaks 
of the Zhang model E { = 4/9 i with i = 0, 1,2 R. 
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FIG. 5. The finite-size scaling plot of the first maximum of 
the energy distribution p(E). The inset displays the statistical 
weights /o and f% as a function of the inverse system size. The 
values of the infinite system (extrapolation to the vertical 
axis) agree with those of the Zhang model (open symbols) 
obtained from If 71]. 



Fig. H we plot p(E) for various values of E c . With increas- 
ing critical energy a peak structure appears with four 
distinct peaks located at Ei = 5/16 i with i — 0,1,2,3. 
Thus the positions of the peaks agree with the corre- 
sponding values of the Zhang model (see JlTj and refer- 
ence therein). 

A similar analysis on a honeycomb lattice is shown in 
Fig. ^. Here the probability distribution is characterized 
by three peaks (z = 3) located at multiples of E = 4/9. 
Again these values are in agreement with those of the 
Zhang model. 

It is known from the Zhang model that the maxima 
of the energy distribution p(E) increase with the sys- 
tem size and that p(E) is characterized by z <5-peaks 
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FIG. 6. The correlation function C(r) for different values of 
the critical energy E c . The inset displays the value C(r = 2) 
which is an indicator for the crossover from the Manna scal- 
ing behavior [C(r — 2) < 0] to the Zhang scaling behavior 
[C(r = 2) > 0]. 




for L — > oo |Tj. Therefore we measured the probabil- 
ity distribution p{E) for a large but fixed critical energy 
(E c = 10 4 ) and for various system sizes. In Fig. || we 
present the rescaled distribution L v p(E) as a function of 
L~ y (E — Ei) around the first peak at E\ = 5/16. The re- 
sulting data collapse shows that similar to the pure Zhang 
model the peak diverges for L — ► oo. But in contrast to 
the pure Zhang model where the finite-size behavior of 
p(E) is characterized by an exponent y w 0.6 (l7]] we get 
for E c — 10 4 the significant different value y = 0.27±0.1. 
Similar results are obtained for the second and third peak 
(not shown) . Further investigations have to show if the 
finite-size scaling exponent y depends on E c , i.e., if y 
tends to 0.6 for E c — > oo. 

But nevertheless the finite-size scaling analysis of the 
probability distribution yields that p(E) is characterized 
by four <5-peaks [Eq. (||)] and one can compare the statis- 
tical weight fi of each peak with the values of the pure 
Zhang model. Therefore we divided the interval [0, E c ] 
in four parts and measured in each part the area under 
the curve p(E) for various system sizes. The values of 
fo and /2 are shown in the inset of Fig. as a function 
of the inverse system size. The same system size depen- 
dence was observed for the pure Zhang model (Tt]]. The 
extrapolation to L — > oo yields the statistical weights fi 
and the obtained values agree with those of the Zhang 
model (see vertical axis of the inset of Fig. ||) . 

Finally we consider in this section the energy corre- 
lations in the steady state. The correlation function is 
defined as 



Civ) 



(E L ' E L '+r) — (Erf) (Erf +r) 

' (E c -l)~ 2 " 



(11) 



In Fig. |(| we plot C(r) along the symmetry axis of the 
square lattice. Starting from the autocorrelation peaks 
at r = the energies of neighboring sites (r — 1) dis- 
play an anti-correlated behavior. This is caused by the 
toppling events and is an typical property of sandpilc 
models where the energy of instable sites is reduced to 
zero. In the inset of Fig. |^ we plot the energy correla- 
tions for the distance r = 2 as a function of E c and K, 
respectively. For small values of E c we get negative cor- 
relations which increases with the critical energy. In the 
surrounding of K* « 0.15 the behavior changes from an 
anti-correlated behavior (C < 0) to a correlated behav- 
ior (C > 0). Positive correlations between sites of the 
same sublattice (r — 2) are a typical feature of sandpilc 
models with an isotropic energy transfer to the nearest 
neighbors. Thus positive correlations within a sublattice 
indicate a Zhang like behavior whereas negative corre- 
lations are a characteristic of the Manna model. The 
analysis of the avalanche distributions in the next sub- 
section reveals that the model exhibits the typical Manna 
scaling behavior (for E c = 2) as long as the energies are 
anti-correlated, i.e., the scaling behavior changes above 
the value K* . 
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B. Avalanche distributions 

In the following we consider the avalanche distributions 
[Eq. (^)] and examine the behavior of the corresponding 
exponents t x as a function of the critical energy. Using 
the reported values of the exponents of the Zhang and 
Manna (for E c = 2) model HJl|Jl|Q we expect that 
the difference is nearly %1 for the size exponent r s , less 
than %2 for r a , and nearly %5 for the radius exponent 
r r . In the case of the duration exponent it seems that 
both models are characterized by the value Tt « 3/2. In 
the following we analysis the size and radius exponents. 
The exponents are obtained from a regression analysis of 
the corresponding probability distribution. 

In Fig. [?] we plot the exponent r r for different values 
of the system size L and different critical energies E c 
(see also Table Q). One can distinguish between three 
regimes: for small values of the critical energy (E c < 
100) the exponents are independent of E c and display no 
significant system size dependence. This is the regime 
of the Manna universality class and the obtained results 
agree with the observed universal behavior for E c < 10 
m|§. 

For large values of the critical energy (E c > 5000) 
another regime occurs where the exponents are nearly 
independent of E c but display a significant system size 
dependence. These system size dependence agrees with 
the observed L-dependence of the Zhang model [Eq. (||)]. 
Thus we recover the third fingerprint of the Zhang model 
for sufficiently large values of E c . We call this regime in 
the following the Zhang regime. Between the Manna and 
the Zhang regime we find a third transient regime where 
the exponents depend strongly on the critical energy and 
on the system size. 

The avalanche size exponent t s displays a similar be- 
havior (see Fig. || and Table |J) ■ Again the Manna regime 
is characterized by independent exponents. Approach- 



K 

1 0.5 0.1 0.05 0.01 o 



1.30 




• L= 128 
■ L= 256 

♦ L=512 
AL=1024 

10 5 ° °= 



FIG. 8. The avalanche exponent r s as a function of the crit- 
ical energy E c and the crossover parameter K, respectively. 
The dashed line indicates the value of the Manna model for 
E c = 2. The corresponding values of the pure Zhang model 
are marked as E c = oo. 



ing a certain value (E c w 100) the exponent depends 
strongly on E c in the transient regime. In the third 
regime (E c > 5000) we observe the characteristic Zhang 
system size dependence. Thus we see that the scaling 
behavior corresponds to that of the Zhang model for 
large but finite values of the energy threshold E c . This 
means that adding stochasticity to the toppling rules of 
the Zhang model does not force a change of the univer- 
sality class in all cases. Only a sufficiently large stochas- 
ticity is relevant, otherwise it is irrelevant and the scaling 
behavior is unchanged. 

Considering the behavior of the exponents r s and r r in 
the transient regime it seems that the exponents increase 
with the system size. Especially the "over-shooting" ef- 
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FIG. 7. The avalanche exponent r r as a function of the crit- 
ical energy E c and the crossover parameter K, respectively. 
The dashed lines indicate the value of the Manna model for 
E c = 2 (upper) and the Zhang (lower) model. The corre- 
sponding value of the pure Zhang model is marked as E c — oo. 



FIG. 9. The avalanche exponent t s as a function of the 
scaling parameter K for the simple cubic (sc) and honey- 
comb (he) lattice. The dashed line indicates the value of the 
Manna model for E c — 2 and the arrow marks the value 
K* — 0.15 obtained from the analysis of the correlation func- 
tion [Eq. (0)]. 
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TABLE I. Some values of the avalanche radius exponent 
r r for various values of the critical energy E c and different 
sizes L of the square lattice. The corresponding value of the 
pure Zhang model is marked as E c — oo. 



_Ec L = 256 L = 512 L = 1024 

2 1.769 1.748 1.737 

5 1.752 1.732 

10 1.747 1.777 

20 1.740 1.735 

50 1.755 1.738 

100 1.752 1.736 

200 1.664 1.651 

500 1.416 1.470 

1000 1.475 1.535 1.581 

2000 1.539 1.624 1.662 

5000 1.551 1.630 1.642 

10000 1.555 1.633 1.680 

50000 1.547 1.634 1.668 

oo 1.539 1.630 1.648 



V. SUMMARY 

We investigated the scaling behavior of the Manna 
model as a function of the threshold condition. Increas- 
ing the value of the critical energy a crossover takes 
place from the universality class where the energy trans- 
fer is undirected on average to the universality class of 
an isotropic energy transfer (Zhang model). For suffi- 
ciently large thresholds all characteristic of the Zhang 
model could be recovered, namely the quantization of the 
energy values in the steady state, the statistical weights 
of the energy quantums in the thermodynamic limit and 
the observed system size dependence of the avalanche ex- 
ponents. Our analysis suggests that the scaling behavior 
of the stochastic sandpile model is connected to the neg- 
ative correlations between next nearest lattice sites. 

I would like to thank A. Hucht for useful discussions. 



feet decreases with increasing L and usually one would 
expect that this regime disappears for L — > oo. Unfor- 
tunately it is impossible to conclude from the numerical 
results whether the transient regime still exist in the ther- 
modynamic limit. Further investigations are needed to 
answer this question seriously. 

In the last figure we plot the size exponent r s for two 
different lattice types as a function of the crossover pa- 
rameter K. As one can see both curves display the same 
behavior, i.e., the transition from the Manna regime to 
the transient regime takes place at K* « 0.2 indepen- 
dent of the lattice type. This value is in agreement with 
the value K* 0.15 obtained from the analysis of the 
correlation function, i.e., the change of the scaling be- 
havior coincides with the change from negative to posi- 
tive correlations C(r = 2). Thus we conclude that the 
Manna scaling behavior is strongly connected to the anti- 
correlations between next nearest neighbors. 

TABLE II. Some values of the avalanche size exponent r s 
for various values of the critical energy E c and different sizes L 
of the square lattice. The corresponding value of the pure 
Zhang model is marked as E c — oo. 



E c 


L = 256 


L = 512 


L = 1 


2 


1.280 


1.278 


1.279 


5 


1.279 


1.276 




10 


1.281 


1.279 




20 


1.282 


1.279 




50 


1.273 


1.276 




100 


1.273 


1.270 




200 


1.224 


1.230 




500 


1.198 


1.198 




1000 


1.219 


1.229 


1.230 


2000 


1.237 


1.242 


1.257 


5000 


1.247 


1.253 


1.261 


10000 


1.246 


1.254 


1.260 


50000 


1.248 


1.255 


1.262 


oo 


1.247 


1.255 


1.260 
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